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ABSTRACT 

Bubbles in the interstellar medium are produced by astrophysical sources, which 
continuously or explosively deposit large amount of energy into the ambient medium. 
These expanding bubbles can drive shocks in front of them, which dynamics is markedly 
different from the widely used Sedov-von Neumann- Taylor blast wave solution. Here we 
present the theory of a bubble-driven shock and show how its properties and evolution 
are determined by the temporal history of the source energy output, generally referred 
to as the source luminosity law, L(t). In particular, we find the analytical solutions 
for a driven shock in two cases: the self-similar scaling L oc (t/t s ) p law (with p and t s 
being constants) and the finite activity time case, L oc (1 — t/t s )~ p . The latter with 
p > describes a finite-time-singular behavior, which is relevant to a wide variety of 
systems with explosive-type energy release. For both luminosity laws, we derived the 
conditions needed for the driven shock to exist and predict the shock observational 
signatures. Our results can be relevant to stellar systems with strong winds, merging 
neutron star/magnetar/black hole systems, and massive stars evolving to supernovae 
explosions. 

Subject headings: ISM: bubbles; shock waves; ISM: jets and outflows 



1. Introduction 

Astrophysical bubbles are formed around sources with outflows in the form of winds, electro- 
magnetic radiation, Poynting flux, etc., when outflows are strong enough to sweep up much material 
in the ambient medium. Examples include supernova remnants, pulsar wind nebulae, stellar-wind- 



driven bubbles around early type stars, Wolf-Rayet stars or star clusters (Gaensler Sz Slane 2006 



Crowthe"r||2007 Maeder Sz Meynet]|2012 ), as well as some models of gamma-ray bursts (Lyutikov & 



^lso at the ITP, NRC "Kurchatov Institute", Moscow 123182, Russia 



Blandford||2003 ) . In a separate paper flMedvedev fc Loeb|2012[ ), we demonstrate that close binaries 



of compact objects (neutron stars, magnetars, black holes) can produce Poynting-flux-driven bub- 



bles during their final inspiral and merger (McWilliams & Levin 2011; Etienne, et al. 2012). The 



bubbles can be accompanied by shocks, which are known to be efficient accelerators of cosmic rays. 
Observationally, the shocks are detected via synchrotron radiation from the accelerated population 
of electrons in the ambient or self-generated magnetic fields, or by Compton up-scattering of lower- 
energy photons by the energetic electrons. Theoretical understanding the dynamics of the bubbles 
and the associated shocks combined with multi-wavelength telescope observations, gravitational 
wave search (e.g., with Advanced LIGO) and neutrino detection experiments (with KamiokaNDE, 
IceCube and others) can provide valuable information about the nature of the central engines. 

In this paper we consider systems in which energy is continuously pumped into the medium 
by a central source. This energy produces excess pressure which pushes on the surrounding plasma 
and results in the formation of an expanding bubble (or a cavity) . If the bubble expansion velocity 
exceeds the sound speed in the ambient medium, then a shock forms ahead of the bubble, as shown 
in Figure [TJ Such a shock wave is different from the Sedov-von Neumann- Taylor blast wave (Sedov 



1946) produced by a point-like instantaneous explosion and freely expanding into the interstellar 
medium (ISM). Instead, the shock is continuously driven by the ever- increasing pressure inside the 
bubble, rendering the classical Sedov-von Neumann- Taylor solution inapplicable. 

Formation and dynamics of astrophysical bubbles has been a subject of intense research (see 
e.g., the review by Ostriker &: McKee||1988 ). However, analytical studies have mostly been limited 



to the self-similar solutions R(t) oc t a , v{t) oc i a_1 with R and v being the size and the expansion 
speed, and a being a constant. This requires the energy deposition rate into the bubble, which 
we colloquially refer to as the "source luminosity" dE^^^/dt = L(t), to be a power-law in time, 
L(t) oc t p , where p is a constant. This is a very restricting assumption because many astrophysical 
sources exhibit an explosive behavior: they are nearly constant-luminosity sources for a long time 
followed by a rapid — "explosive" — increase in energy output, formally tending to infinity within 
a finite timed 

L(t)=L s (l-t/t s )-P, (1) 

where t s is the time during which the source is active, i.e., its "lifetime", L s is the source luminosity 
at early times t <C t s and p is the luminosity index, determining the rate of energy deposition in 
the bubble. Such a luminosity law with the finite-time singularity (FTS) is not rare in nature. 
For example, a binary consisting of two compacts objects, i.e., a neutron star (NS), a magnetar, 
or a black hole (least one of the companions should be magnetized) is producing a bubble filled 
with material with the relativistic equation of state ( Medvedev &: Loeb]|2012 ) and its luminosity is 



described by Eq. ([!]) with the index p being around 3/2; the exact value depends on specific details 
of the energy extraction. The merger of two neutron stars is considered to be the conventional model 



1 The singular (explosive) behavior occurs if p > 0. However, our analysis below makes no assumption about the 
value of p. 
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of short gamma-ray bursts (GRBs). However, the formation of the bubble and the associated driven 
shock and their observational signatures have not been addressed. The bubble-driven shock can, 
however, produce a unique precursor to the main short GRB event and serve as a benchmark of the 
NS-NS progenitor. Combining these results with the detection of gravitational waves accompanied 
by a neutrino signal measured, respectively, by Advanced LIGO and other gravitational wave 
detectors (during the last moments of the binary inspiral) and by KamiokaNDE, IceCube and other 
neutrino detectors (during the merger event itself) will also greatly advance our understanding. 
Thus, the development of the theory of the bubble and shock dynamics in such systems is of great 
importance for astrophysics. 

The paper is organized as follows. In Section [2] we develop the theoretical model of a bubble 
and the bubble-driven shock and obtain the master equation describing the system evolution for 
an arbitrary luminosity law L(t). In Section [3] we find analytical solutions for the master equation 
for both the self-similar and singular (i.e., FTS) luminosity laws and explore the conditions under 
which the driven shock solutions exist. We also determined the conditions under which the exact 
FTS solutions (for the radius and velocity) exist. However, realistic systems, such as the neutron 
star binaries, do not fall into that category. Nevertheless, we were able to find simple analytical 
approximate FTS solutions for them as well. Finally, we demonstrate that the obtained analytical 
solutions agree very well with the full numerical solutions in the appropriate limits. Finally, in 
Section [4] we present conclusions. 



2. The bubble+shock model 

The system at hand is depicted in Fig. [T] For simplicity, we consider spherically symmetric 
systems. The bubble is filled with some material with pressure pbubbie arid the equation of state 
parameterized by the adiabatic index %, which depends on the bubble nature and composition. If 
the bubble is driven by Poynting flux (e.g., in a form of a strongly magnetized relativistic e~e + 
wind) as in pulsar wind nebulae, or double neutron star or magnetar binaries, then the value of the 
index is % = 4/3. If the bubble is driven by a non-relativistic strong stellar wind, then jf, = 5/3. 

The bubble surface acts as a piston and pushes the ambient plasma outward to produce an 
outgoing shock wave. Is the outgoing shock relativistic? Let us assume, for example, that the total 
energy E ~ 10 46 erg is deposited into the bubble over some time t souvce during which the source is 
active. If the shock is relativistic, it should expand with the speed ~ c and the Lorentz factor 

EE 

r ~ — ~ » i (2) 

Mc 2 (47r/3)(ni SM m p )(ci 

source J <-- 

where M is the mass of the ISM material swept by the shock over time £ S ource and m p is the proton 
mass (assuming purely hydrogen ISM). For the ISM with a uniform density nisM ~ 1 cm -3 , we 
have the constraint on the source lifetime: 

Wee «4x 10 5 £^ 3 n^ 3 s ~ 4.5 n KMo da y s > ( 3 ) 
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where E^ = E/(10 AG erg), nisM,o = ^ISm/(1 cm" 3 ) and similarly for other quantities henceforth. 
Thus, if the typical source activity is much longer than a few days, the shock can be assumed to 
be non-relativistic. This constraint is even less restrictive for lower-energy sources. 

The shock is propagating into the unperturbed unmagnetized or weakly magnetized ISM with 
mass density pism and the adiabatic index 7 = 5/3. For simplicity, we assume the ISM to be cold, 
so that we can neglect its pressure; hence the shock is a high-Mach-number strong shock. Besides, 
weak shocks are less interesting from the observational point of view anyway: being just a mild 
perturbation to the medium, their observational signatures are hard to detect. Finally, between the 
shock and the bubble lies the shell of the shocked gas, which mass density and pressure, /9 s heii an d 
Pshelh are determined by the Rankine-Hugoniot shock jump conditions. The pressure equilibration 
time behind the shock is assumed to be fast enough to establish pressure equilibrium throughout the 
system, hence the bubble-shell interface plays a role of a contact discontinuity and ptmbbie = Pshell- 

We assume the central engine deposits energy (e.g., electromagnetic, kinetic, etc.) in the 
bubble with luminosity Lit) = dE/dt, which is a function of time. This power goes into the 
internal energies (i) of the bubble, dUbubbie, and (ii) of the shocked gas shell, dU s h e \i, (hi) the 
change of the kinetic energy of the bulk motion of the shell, dK s ^ e n, assuming its swept-up mass 
Mswept = const, (vi) the change of the kinetic energy, <ifr@shock; of the newly swept gas, dM swept , 




Fig. 1. — Schematic representation of the system. The central source outflow produces a bubble, 
which creates a shock in the ambient medium and the shell of the shock-heated gas. 
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and also (v) heating of this shocked gas, cii7@ s hock> to satisfy the Rankine-Hugoniot relations. We 
can neglect the pdV work due to the expansion because the external pressure is assumed to be 
vanishingly small in the cold ISM. Thus, the Master equation is 

L(t) dt = (it/bubble + dUsheW + dK^w | M t + dU@ s hock + d-K@shock- (4) 

In addition to this equation, the Rankine-Hugoniot relations for a strong shock should be used. 
The continuity equation yields 

Pshell _ Ui _ ^shock 7 + 1 _ 

— — 7~B \ — 7" = K W 

PISM U2 ^shock - 'Wshelll-Kshock) 7 - 1 

and 

^shell-shock) -. -1 , e s 

= 1 - « , (6) 

^shock 

where k is the constant compression ratio of a strong nonradiative shock, u\ and U2 are the upstream 
and post-shock velocities in the shock comoving frame which transform to the lab frame as u s hock = 
u\ and i>shell(-Rshock) = U1-U2. From the momentum conservation, Pism + Pism^i = Pshell + Psheliuj , 
for a cold ISM (i.e., neglecting pism), one gets 

Psheii = PiSMt4 ock (! - • ( 7 ) 
We now calculate the terms in Eq. Q. 



2.1. Internal energies 

The internal energies of the bubble and the post-shock gas shell are those of ideal gas: 

^bubble = d( Pbubble^bubble ) , (8) 

V76 - 1 / 



dU s hel\ = d 



r PshellKhell J , (9) 

7-1 J 



where ^bubble = (47r/3) -Rabble and ^eii = Khock - ^bubble are the volumes of the bubble and the 
shell, respectively, and V s hock = (47r/3)i?g hock and recall that ^bubble = Psheil- The swept- up ISM 
mass occupies the post-shock shell, M swept = /OisM^hock = PshellKshell; hence one has the relations: 

^shell -1 ^bubble , -1 , ln N 
= K , — = 1 - K . (10) 

^shock ^shock 

Therefore, 

^1 t\j ^ / 2 \ 

"^bubble = — _ ^ PISM d (Wshock^shock) , (H) 

dU shc \\ = _ j- PISM d (^shockKihock) • (12) 
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2.2. Kinetic energy 



To calculate the shell kinetic energy one needs to know the post-shock velocity profile v s heii(-R)- 
The post-shock gas is hot, hence thermal conduction is rapid and the temperature can be taken 
uniform throughout the shell. Since both temperature and pressure are uniform, then density is 
uniform as well, p(R) = p s hell = Pismk = const. Continuity equation of an incompressible gas in 
a steady state, V • v = 8r (R 2 v s i ie n(R)) = 0, together with the boundary condition v (i? s hock) = 
Wsheii(-Rshock) = ^shock (1 - k™ 1 ) yields 

/ hi ; i \ 

— ^shell (-^shock) 



(-^shock 



(13) 



The kinetic energy of the shocked gas (see Fig. [T]) is 

M swept dMVshcn ( R )2 



K 



shell 



-^swept^shock 



(14) 



where 



3r (1 

12 

f 7 " 



-i\ 2 



7 + 1 



1 - K 

1/3 



-1/3 



6 2/3 



and we have used that dM = (/jism«) <±irR 2 dR and i? b ubbie/ ^shock = (^bubble/Kshoek) 
At last, the change of the shell kinetic energy due to acceleration/deceleration is 



1/3 



G^shelll 



swept 



PiSM(4vr/3)^ hock C d(v 2 shock /2) . 



(15) 
(16) 



2.3. Shock contributions 

As the shock propagates, it heats up the cold ISM gas of mass dM swept = pism ^Kihock to the 
post-shock temperature T s h e ii = Psheii/^sheii = m p u shock ~~ ft -1 ) k ~ 1 ano - a ^ so accelerates it to the 
post-shock velocity, i> s heii(^shock) = I'shock (l — k -1 ). Thus, at the shock 

(l-K- 1 )^ 1 2 
d^Whock = ; PISM ^shock^shock' (l?) 

7-1 

iP 1 2 

d-K@shock = PlSM^Kshock S °° (l — «T ) . (18) 



2.4. Final analysis 



Equation Q together with equations (11), (12), (16), (17) and ( [181 ) determine the evolution of 
a shock driven by the pressure inside a cavity. Given the luminosity, L(t), it allows us to determine 
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-Rshock(^)) because the shock velocity is w s hock = dR^o^/dt = -R s hock and the shocked volume is 
Kihock = (^ 7r /3)^shock- ^11 other quantities, e.g., -Rbubbie, Pshellj etc. follow straightforwardly from 
the equations above. Hereafter we will often omit the subscript "shock" wherever it does not cause 
a confusion. Upon the substitutions, Eq. Q becomes 



L{t) 



(4/3)7rp IS M 



d(R 3 R 2 ) t^dR 2 .- 2 dR 3 

1 ^^ + 2 R ^r +CR ~& 



(19) 



where 



(1 



.-l 



76-1 
2(7fe + 1) 
(7 + l) 2 (76- 



+ 



7-1 

63 
~ 32' 



(20) 



1 



7 



12 



1 

27 
16' 



+ 



- ( 7 + l) 2 ~16- (21) 
Here the first term represents the internal energies of the bubble and shell, the second term is the 
bulk kinetic energy of the shell and the last term is the contribution of the shock. Hereafter we 
use, for concreteness, that the bubble is filled with relativistic material (e.g., a relativistic plasma, 
or a highly magnetized wind, or electromagnetic radiation) with 75 = 4/3, the ambient gas is 
non-relativistic with 7 = 5/3 and the compression ratio is k ~ 4. 



Master equation ( 19 ) is the inhomogeneous second-order nonlinear differential equation, which 



(3?/ + C)R 2 R 3 + (277 + £)R 6 RR. 



(22) 



can further be simplified to yield 

Lit) 
(4vr/3)pi SM 

This is the main equation of our analysis. For any function of the source emission luminosity, 
L(t), pumping energy into the bubble, this equation describes the evolution of the shock radius, 
-^shock = R(t), and the associated parameters, e.g., the shock velocity, v(t) = R(t), the size of 
the bubble, i^bubble(i) = (l — k -1 ) 1 ^ 3 R{t) ~ (3/4) 1 / 3 ii(t), etc. We stress that this equation is 
applicable to a bubble of any origin. For example, it describes a Poynting-flux-driven bubble 
formed by an inspiraling binary; in this case 7^ = 4/3. It also describes a bubble blown by a strong 
stellar wind with the kinetic energy luminosity given by L{t) and the adiabatic index of the gas in 
the bubble 7^ = 5/3. 



3. Results and observational predictions 

To proceed further, one needs to specify the source luminosity law, L(t). Below, we will consider 
the two most common ones: the self-similar law L(t) oc t p and the "explosive" finite-time-singular 
law L(t) oc (t s - t)~ p with p > 0, 
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3.1. Self-similar solution 



3.1.1. Structure 



We first look for a self-similar solution. Let us assume that the source luminosity and the 
shock position are power-law functions of time 



L(t) = L s (t/t s 



R(t) = R s {t/t s 



(23) 



where R s and a are constants to be determined and L s , t s and p are known constants set by the 
source physics. This self-similar solution is valid for the duration of the source activity, i.e., while 
the bubble pressure is high enough to push the shock; at (much) later times, the shock dynamics 
should asymptote the Sedov-von Neumann- Taylor solution. 



The self-similar solution (23) to equation (22) is 

P + 3 



a 



L s t s 



1/5 



4vr pisM-V 



(24) 



where 



A 



(3r ? + C)a 3 + (2r ? + ^)a 2 (a-l) 

5 6 2/3 + ff + 6 2/3 

4 V 8 



9 



-a 



(25) 



The scalings given in equations (23), (24) agree with those in ( Ostriker McKee||1988 ). The above 
solution is meaningful if R > 0, i.e., A > and, hence, 



« > "crit 



6 2 / 3 - 5/4 
6 2 /3 + 17/8 



and 



V > Pcrit = 5a c 



0.378 



■1.11. 



(26) 



(27) 



This condition means that the energy injection in the system cannot be too slow; otherwise the 
shock would move too fast for the contact discontinuity to catch up with it and the assumption of 
the pressure equilibration breaks down. 



3.1.2. Estimates 

We have obtained that given the luminosity of the source, L(t) = L s (t/t s ) p , the shock evolution 
is given by R s hock(t) = R s (t/t s ) a and v shock (t) = (aR s /t s )(t/t s ) a ~ 1 with a = (p + 3)/5. The value 
of R s is rather insensitive to the numerical value of A (unless a is very close to the critical value), so, 
R s ~ (L s tg/ pism) 1//5 is a rather good order-of-magnitude estimate. More accurately, assuming the 
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source activity to be of the order of a hundred days, t s ~ 10 s, and the luminosity L s ~ 10 erg s , 
we estimate that 

Rshock(t s ) = R s ^ 4.3 x 10 16 A-^Ll^t^n'^ cm, (28) 

v sho ck(t s ) = aR s /t s ~ 4.3 x 10 9 ai^/^^^fn^ cm s" 1 . (29) 

If a < 1, then the shock velocity oc is decreasing with time, thus at some earlier time it was 
~ c, so the non-relativistic approximation was invalid. For the solution to be valid at early times 
for the duration of the source activity, the shock velocity should be increasing with time, i.e., a 
should be greater than unity, hence the luminosity index must be p > 2. 



3.1.3. Observable signature 

We assume that the shock accelerates electrons via Fermi process and generates/amplifies 
magnetic field. The relativistic electrons in magnetic fields produce synchrotron radiation which 
can be observed. We assume that the electrons and magnetic fields carry, respectively, fractions e e 
and 6b of the internal energy density of the shocked gas, cf., Eq. (12), 

Mshcll = 7- PlSM^ shock , (30) 

7-1 

that is 

7 e m e c 2 n eiShe n = e e u s heii, B 2 /8tt = e B u s h e ii, (31) 

where j e m e c 2 is the average energy of an electron and n e)S h e ii = nne,iSM — k^ism * s the number 
density of electrons in the shocked gas shell. Numerically, the average Lorentz factor of accelerated 
electrons is 

_ ,.x k-1 m p fv shock \2 

7e(t) = ee ^(^T)m e {— ) 

* 11 Uc?A-^L%t^ri^^\ (32) 

so the bulk electrons are mildly relativistic 7= ~ 3 for a typical acceleration efficiency of e e = 0.3. 
If the radiating electrons are distributed in energy as a power-law with index s as dn e jd^ oc 7~ s 
with a minimum Lorentz factor 7 m , then 7 m = 7 e (s — 2)/(s — 1). The magnetic field strength is 



B (t) = ( g B ^~ _^ 87rm p rai SM t4 ock ^ 



1/2 



x 3.0 x 10- 2 efaA-y%%t:fn%l Q tr l G, (33) 

that is, B is of the order of a milligauss for a nominal 6b ~ 10~ 3 , which is much larger than the 
typical microgauss fields in the ISM. This means that the field should be generated or amplified at 
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the shock by an instability (Weibel, Bell, cyclotron, and such), or by MHD turbulence, or by some 
other mechanism. 

Relativistic electrons emit synchrotron radiation. The peak frequency of the radiation spec- 
trum is in the radio band, 

m e c 

~ 9.3 x 10 6 &J*cl*A-\8 - 2)/(s - l)] 2 X Sj39 t-?nr s i;>? (a - 1) Hz - (34) 

The total emitted power by a relativistic electron is P = (4/3)<7tC7 2 (B 2 /8ir), where ot is the 
Thompson cross-section. Thus the spectral power at the peak (measured in erg s^ 1 Hz -1 ) is 
Pv,ma.x(t) ~ P/v s — (crTm e c 2 / 3e) B . The observed spectral peak flux from a source located in 
our galaxy at the distance D = 10 23 cm (i.e., ~ 30 kpc) is F Vjmax = N e P U;max /(4-KD 2 ), where 
N e = n. e> isM^shock is the total number of emitting electrons, hence 

1 f 4tt 3 \ aTm e c 2 

Fu,max{i) = [ "y^shock^ISM ) ^ B 

* 3.0 x 10 3 D^c^L^r^T 1 (35) 

Although the peak frequency can fall in the self-absorbed part of the spectrum, as is typical of su- 
pernova shocks (e.g., Waxman Sz Loeb|1999 ), the peak flux above is still useful for the normalization 
of the spectrum 

/ v \-(-l)/2 

F v {v,t) = F u>max (t) i^j^J 

K ^-( S -l)/2 ^ (3-5 S +a(3+5 S ))/2 

OC ^-0-75^7.750-4.75 ^ (36) 

the latter scaling corresponds to the nominal value of s = 2.5. We remind that the index a is 
related to that of the energy injection L(t) oc t p as a = (p + 3)/2. 

Let us write F v (u, t) oc u a t b , where a and b are spectral and temporal indexes determined from 
observations. Then one can infer the (effective) energy injection index as 

p = (l-5o + 6)/(4-5o). (37) 

Also, for a non-self-similar luminosity law, one can define the effective index as 

p cS = dlogL(t)/dlogt, (38) 

which can be compared with observations. 

Finally, we estimate the self-absorption frequency as the frequency at which the optical depth of 
the emitting region of thickness -R s hock — -^bubble is of the order of unity, t u = k v (i? s hock — ^bubble) ~ 1 5 
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where k v is the self- absorption coefficient (Rybicki & Lightman 1979): 



8 ™^7 whh ihl?wJ- (39) 

A power-law distributed electrons dn e /d^ = (n eiS h e ii/7m)(7/7m)~ s produce the spectrum (per elec- 



tron) Pud) — Pi/^axfiv/vc), where f(x) is the function F(x) of Rybicki Sz Lightman (1979) 
involving the integral of a modified Bessel function, up to a numerical factor. Using a dimensional 
analysis, noting that n ejS h e ii(-Rshock — -^bubble) — ^iSM-Rshock/3 and absorbing all numerical factors 
into a constant K ~ 10 -2 , one gets the absorption frequency, v a , 

v a ~ (Ka T c%n 1S MRs^/me) 2/{s+4) v^'^. (40) 

which depends very weakly of the constant K, the density and the shock radius and is nearly 
independent of the synchrotron frequency for typical spectral indices s > 2. For typical values 
of the parameters, v a is of the order of ~ 10 8 — 10 9 Hz. Note that this analysis assumes that 
all the electrons are relativistic; otherwise, if only a fraction r] ie \ < 1 of them is relativistic, the 
self-absorption frequency is lower by ~ 7 lrl\ S+4 ''' '■ 



3.2. Solution with a finite-time singularity 

3.2.1. General structure 

Traditionally in astrophysics, one looks for a self-similar solution, which is presented above. 
However, for some astrophysical systems, self-similarity may not be a good approximation. For 



example, in the paper (Medvedev Sz Loeb 2012) we have shown that the evolution of the binary 



system is described by a finite-time singular (FTS) solution, because the inspiral takes a finite time. 
We have found that the luminosity is 

L(t) = L s (l-t/t s )-P% (41) 

where L s , t s and p* are known "source" constants set by its physics and initial conditions. To be 
precise, these scaling should break down at some earlier time t m < t s , otherwise the total energy 
diverges if > 1, which is the case of a neutron star binary, for example. It also breaks down 
on physics grounds because general relativistic and other effects, as well as finite object sizes, were 
omitted from consideration. The above form of L(t) suggests to look for a solution in the form 

R(t) = R s 4l-t/t s )~ a % (42) 

where R s and a* are constants to be determined. Such a solution is formally valid at times t < t s . 
At late times t 3> t s , there is no energy injection, therefore the shock dynamics should asymptote 
the Sedov-von Neumann- Taylor solution. 
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The finite-time singular solution, Eq. (42), to the Master equation (22) is 

3 L s £ ^ 1/5 



a* 



p* - 3 



Rs,* 



4vr pism A, 



(43) 



where 



(3T7 + C)a2 + (277 + Oa2(a. + l) 



-^ + 6 2 / 3 + a^ + f> 2 '* 



The solution is physically meaningful if R > 0, hence A* > 0, 

5/4 -6 2 / 3 



Ck* ^> a* crit 



and 



17/8 + 6 2 /3 



0.378 



(44) 



(45) 



P* > P*,crit = 5Q* >crit + 3 ~ -1.11. (46) 

This condition constrains the energy injection index: if < p* jCr it then the energy injection rate 
is not enough for the contact discontinuity to catch up with the shock. 

Unlike the self-similar solution, this solution with a finite time singularity has another con- 
straint. The physically meaningful solution (in this particular class of solutions) is the one of the 
expanding shock; hence 

a* > and > 3. (47) 

We emphasize that the above scalings are applicable for the duration of the central engine activity 
only, t < t s . Moreover, the solution discussed in the previous subsection does not have much 
physical meaning, because the system size increases to infinity in a finite time. Perhaps, it can 
make sense at times substantially smaller than t s , when relativistic and other effects are negligible. 

As we mentioned earlier, some astrophysical systems, such as the neutron star and magnetar 
binaries (Medvedev & Loeb 2012), have the luminosity indexes p* = 3/2 or 7/4, depending on the 



electromagnetic energy extraction mechanism. For suchp* values from in the range, p* :C vit < P* < 3, 



the bubble-driven shock exists, but the solution given by equation (42) describes an unphysical 
collapsing shock with R(t) — > as t — >• t s . This result is independent of most of the assumptions 
and easily follows from the energetics considerations or just the dimensional analysis, cf. equation 



(22) 



(1 - t/t s )- p * oc L(t) ~ d{Mv l )/dt ~ pd(R 3 R 2 )/dt oc (1 - t/t s ) 



-5a, -3 



(48) 



which yields the relation = 5a* + 3. On physics grounds, for values p* iCr it < P* < 3, the 
shock radius tends to a constant, R — > i? max as t — > t s , but the form (1 — t/t s )~ a * does not 
have such asymptotic for any a* ^ 0. A different form is needed, but it seems unlikely that one 



can find a general or partial solutions to Eq. (22) directly, because the solution should be of the 
form J dt\J (1 — t)~P* +l + C, which cannot be expressed in elementary functions for an arbitrary 
p*, though it can be expressed in elliptic functions for some rational p*. Below we obtain an 
approximate solution for such a regime. 
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3.2.2. Approximate realistic FTS solution 

In the previous section, we have found that the bubble-driven shock can exist in systems with 
the luminosity index > p* )Cr it! but if is in the range p* jCr it < p* < 3, the physically meaningful 
solutions are not described by the pure FTS solutions. Instead, the physical FTS solution should 
have the property that -R s hock — > -Rmax = const, but v s hock oc (1 — t/t s )~ a ~ 1 — > oo as t — > t s , with 
a* being a new constant. The general and/or exact analytical solution of this kind does not exist, 
so we construct here a composite approximate solution. 

First, we notice that the dependence L(t) = L s (l — t/t s )~ p * implies a constant luminosity 
source at early times, t <C t s , that is L(t) ~ L s (t/t s )°. This case corresponds to the self-similar 
solution with p = in Section 3.1 for which the solution exists: 

# s hock(i) = Rs(t/t s ) 3/5 -> R s , as t-H„ (49) 



where R s is given by equations (24) or ( |54[ ). Second, we make substitutions 

R = R S , R = v shock (t)=v s (l-t/t s )- a *-\ (50) 



where v s is a constant to be determined. By examining of the right hand side of equation (22), 
one sees that both terms are divergent but the last term dominates as t — > t s , because < 3 and 
hence a* < 0. Keeping the leading term, one has 

* _ P* ~ 3 /3 L s t s Y /2 = ( A V /2 jg 



1/2 

where A ~ 4.3 is given by equation (25) with a = 3/5 and the prefactor [A/(2n + £)] — 0.97. 

Thus, we have found an approximate solution describing evolution of the bubble-driven shock 
at the late times: 

v shock (At) ~ v s (At/t s )-^-^ 2 , (52) 
^shock(At) ~ R s (At/t s )° = const, (53) 

where we introduced a new variable At = (t s — t)—>0, which is more convenient when t < t s . Note 
that the obtained regime is rather interesting: the shock is rapidly accelerating but its size and 
swept-up mass remain nearly constant. 



3.2.3. Estimates 



The characteristic values of the shock radius and velocity in the solution above are given by 
equations (52) and (53). However, these values follow from an approximate analytical analysis. 
Comparison with the exact numerical solutions (Section 3.3) yields ad hoc correction factors for R s 
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and v s , see Eqs. (65) and (66), namely xr — 1-33, Xv — 0.51. We use these values in the estimates 
below. We have 

R s * 3.2xl0 16 x/?OKVsM 5 oCm, (54) 
v s ~ 3.1 x 10 9 Xv li^tjj /5 n^ cm s -1 , (55) 

where we assumed 'jb = 4/3 and 7 = 5/3, for concreteness. We, thus, have 

v shock (At) ~ 3.1 x 10 8+ ^ X^ s V 39^1 9 " 5P * )/10 Vsm 5 oA^ (P ^ 1)/2 cm s" 1 

~ 9.8 x 10 9 xv Ll%t: 3 7 /2 \^ At~^ cm s" 1 , (56) 

where we used a nominal value p* ~ 3/2. Note that if p* > 1 then t> s hock oc {At/t s )~^ p *~ 1 ^ 2 — > 00, 
as At —7- and w s hock approaches the speed of light at times 



At < 10 7-2/(p.-l) t ( L l/5 -2/5 -l/5\ 2 /(f* !) 



(57) 



that is, about 10 3 s before the "explosion" time t s , for p* = 3/2. At these times our assumption of 
the nonrelativistic shock breaks down and a different analysis is needed. Note also that at this time, 
the dynamical time of the bubble needed to establish pressure equilibrium throughout is longer, 
R s /c ~ 10 6 s, so an accurate analysis should instead be involved to find a full dynamical solution. 
Such a consideration goes beyond the scope of the present paper. 

Regardless of the model assumption used, plasma processes and details of particle acceleration 
impose additional constraints as follows. If the characteristic dynamical time of the system, which is 
~ At, is longer than the inverse collisional frequency u~r, then a collisional shock forms. Otherwise, 
when At < v~y, then the shock is collisionless. In this case, if At is (much) shorter then the 
Larmor frequency in the ambient field, the shock structure is not sensitive to the ISM field, but, 
instead, is determined by kinetic plasma instabilities (e.g., electrostatic Buneman or electromagnetic 
Weibel-like ones driven by particle anisotropy at the shock). The shortest associated timescale is 
the ion plasma time ~ 10 3 n.^/ 2 s and, moreover, it takes about a hundred seconds for 
an electrostatic shock [or (Vshock/ c ) seconds for a Weibel shock] to form and respond to the 
changing conditions; it takes even longer for particle Fermi acceleration. Thus, the synchrotron 
shock model should be used with great caution for At as short as tens of milliseconds or less. 



3.2.4. Observables 



The above scalings can be used in equations for observables in Section 3.1.3 (note, most of the 
parameters turn out to be functions of the shock velocity alone). In particular, for At ~ 10 5 s, i.e., 
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about ~ 1 day before a merger of an explosion, we have 

>v (At) ~ 5 5 X 1(r 2 + 2 P' v 2 f T 2 /5 f -9/5+ P , -2/5 . -(p.-l) 

~ 55 v 2 e L 2/5 r 3/1 V 2/5 At" 1/a f 58) 

B(At) * 2.1 x 10-^ X ,e^^ )/M ^oAt^- 1)/2 G 

* 6.7 x 10- 4 /2 ^ 3 k1 /2 ^^5 V4 G > ( 59 ) 
*(Af) * 1.8 x 10^ xl 4e¥l(s ~ 2)/(- - l^M^^^n^AtB Hz 

* 5.7 x 10 8 X * ele^Ks - 2)/(s - l)] 2 L Si39 t^\-^ Atf /A Hz, (60) 

P fA^I ~ gn Y ml+l'. v 3 v r)-2 1/2 j-4/5 .(9+5p»)/10 7/10 (p«— 1)/2 j 

^i/,max(At) ~ 8.9X10 XrXv^ 2 3 € B L s,39 t s,7 "lSM,0 A *5 J y 

~ 2 8 x 10 3 y 3 v D- 2 e 1/2 L 4/5 i 33 / 20 n 7/10 At~ 1/A Tv (61) 

F„(At) OC I/ -(-l)/2 At -(5»-3)(p.-l)/4 

OC I/ -(-l)/2 Ar (5--3)/8 < (62) 

Here we remind that one should not use these scalings for too short At, as is discussed in the end 
of the previous subsection. 

One can reverse the argument here and ask: What physical parameters of the system can be 
inferred from observations? Obviously, if the spectral slope and the light-curve indexes of the flux 

F v (t) <xiA(At) A (63) 

are measured in observations, one can readily determine the energy injection index p*: 

1 - 5/3, - 2A 

P * = 1-5/3, ' (64) 
where f3 u ^ 1/5, otherwise s = 3/5 and fit = 0. 



3.3. Comparison with full numerical solution 



The full numerical solution of equation ([22]) and the analytical solutions described in previous 
sections are shown in Figures [2] and [3j Figure [2] shows the shock radius as a function of time. 
The solid curve is the exact numerical solution with p = 7/4 and the dashed curve represents the 
approximate self-similar solution given by equations (23), (24) with index p = 0, which corresponds 



to the early-time asymptotic of the realistic evolution (see discussion in Section 3.2.2). The agree- 
ment of the self-similar solution with the exact one is remarkable. The noticeable deviation occurs 
only at the very late time, just before the coalescence time, but even then the difference is within 
a factor of order unity. Thus, the assumption of i? s hock — > const, as t — > t s , used in Section [3.2.2 



to derive the FTS solution, is justified. We determine numerically that the analytical solution 
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underestimates the final radius (at t = t s ) of the shock by a factor of 

-Rs,exact/-Rs,selfsim = Xi? ~ 1-33 (65) 

which we have included in the estimates of shock parameters and observables. 

Figure [3] shows the shock velocity as a function of time: t/t s (left panel) and of the time before 
singularity At/t s = (t s — t)/t s (right panel). The former elucidates the agreement of the exact 
numerical solution with the self-similar solution and the latter - with the FTS analytical solution. 
Here the solid curve is the exact numerical solution, the dashed curve represents the approximate 
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Fig. 2. — Log- log plot of the shock radius driven by an expanding bubble as a function of time. The 



solid curve represents the full numerical solution of Eq. ( 22 ) and the dashed curve is the self-similar 



solution given by Eqs. (23), (24) 
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Fig. 3. — Log-log plot of the driven-shock velocity as a function of time: t/t s (left panel) and 
1 — t/tg (right panel); note that in the right panel, time evolution is from-right-to-left: the merger 
occurs at At/t s — > 0. The solid curve represents the full numerical solution of Eq. (22), the dashed 



curve is the self-similar solution given by Eqs. (23), (24) and i> s hock = R(t), and the dot-dashed 



curve is the approximate finite-time-singular solution given by Eqs. (51), (52) with a correction 



factor included, see Eq. (66) 
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self-similar solution given by equations (23), (24) with index p = and the dot-dashed curve 



represents the approximate FTS solution. This analytical solution given by equations (51), (52) is 
found to overestimate the velocity by a factor of 2, i.e., 



fs,exact/^s,FTS = Xv 



0.51. 



(66) 



Thus, the analytical solution represented in this figure is the one given by Eq. (52 ) with v s replaced 



with XvV s - Note the remarkable agreement between the exact numerical and analytical solutions. 



4. Conclusions 

In this paper we considered the formation and evolution of an astrophysical bubble and the 
bubble-driven shock propagating in an ambient ISM of uniform density under the assumptions of 
spherical symmetry of the pressure balance throughout the system. The equation of the dynamics of 
the shock (and all other parameters) has been accurately derived for an arbitrary energy rate output 
by the central source, which we colloquially refer to as the "luminosity law," L(t). Furthermore, 
we derived the analytical solutions for two special cases: the self-similar scaling, L(t) oc t p , and 
the finite-time-singular case, L(t) oc (t s — t)~ p , where t s is the source life-time, and p is a constant 
index. The latter "explosive" law can represent, for example, the energy output from a merging 
neutron star or magnetar binary. The analytical solution of the finite-time-singular type has not 
been derived before. 

We have found that the dynamics of the bubble-driven shock is markedly different from the 
classical Sedov-von Neumann- Taylor solution even in the case of the "explosive" FTS behavior. 
The driven shock solutions only exist if the energy injection rate is not too low, namely the index p 
shall exceed some critical value, p > p c , which depends on the equations of state of the bubble and 
the ISM. In particular, for the standard ISM with the adiabatic index 7 = 5/3 and the bubble filled 
with the material with the relativistic equation of state (magnetized plasmas, electron-positron 
plasmas, electromagnetic radiation), % = 4/3, the critical value of p c is —1.11, for both types of 
solutions. Otherwise, if p < p c , the bubble expansion is too slow to catch up with the outgoing 
shock. For p > p c , the self-similar and FTS solutions for the shock radius and velocity are 

^scifsim(i) oc &+ 3 ^ 5 , v sclisim (t) oc t (p " 2)/5 . (67) 
R FTS (t) oc (t s - t)^/\ v FTS (t) oc (t s - t)(^ 8 )/ 5 . (68) 

This FTS solution is physical only if p > 3; otherwise it is unphysical because it describes a 
converging shock. For — 1.11 < p < 3 we have found a physically meaningful approximate solution 

«FTs(i) oc {t 8 - t)°, v FTS (t) oc {t s - t)-^/\ (69) 

Interestingly, the derived solutions are also applicable to the class of systems with finite-time 
but non-singular luminosity laws: L{t) oc (t s — t)~ p with —1.11 < p < 0. These are the systems 



- 18 - 



which have a declining with time energy deposition rate and which are active for a finite time t s . 
In the systems with p > 1, the shock is accelerating as t — > t s . Therefore, the bubble-shell interface 
may become Rayleigh- Taylor unstable if a less dense plasma of the bubble is pushing on denser 
ambient medium. Strong mixing is expected in this case. The overall dynamics may somewhat be 
affected, but the salient features should remain the same because the assumption of the pressure 
balance still holds. Another limitation of our analysis is the neglect of relativistic effects, which 
may be important if the shock velocity is singular at t — )■ t s . We have also used the strong shock 
approximation throughout the analysis. 

Finally, we made observational predictions. We calculated the emission light-curves of the 
bubble+shock systems for both self-similar and FTS cases. We predicted that one can deduce the 



luminosity law index p from the spectral and temporal indexes, see Eqs. (63), (64). Our results 
may be relevant to stellar systems with strong winds, merging neutron star/magnetar/black hole 
systems, as well as massive stars evolving to supernovae explosions. 
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